Virial expansion coefficients in the harmonic approximation 
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The virial expansion method is applied within a harmonic approximation to an interacting N- 
body system of identical fermions. We compute the canonical partition functions for two and three 
particles to get the two lowest orders in the expansion. The energy spectrum is carefully interpolated 
to reproduce ground state properties at low temperature and the non-interacting large temperature 
limit of constant virial coefficients. This resembles the smearing of shell effects in finite systems with 
increasing temperature. Numerical results are discussed for the second and third virial coefficients 
as function of dimension, temperature, interaction, and the transition temperature between low and 
high energy limits. 
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I. INTRODUCTION 

The virial expansion is a classical concept [IHl] which 
has been extended to be applicable for quantum mechan- 
ical systems Q. The expansion is in terms of few-body 
correlations and therefore most efficient when the influ- 
ence of TV-body effects decrease with N. In practice, this 
decrease has to be very fast, because higher order correla- 
tions are extremely difficult to obtain by accurate calcu- 
lations. This fact is not obvious since only the spectrum 
of ./V interacting particles is needed, not wave functions or 
structure nor any other properties. However, obtaining 
these spectra imply solving the N-hody problem which 
is already demanding beyond two particles for general 
interactions. 

In the classical textbook by Huang on statistical me- 
chanics [f| , the virial expansion is discussed for the quan- 
tum mechanical case and it is elegantly demonstrated 
how the second virial coefficient can be obtained from 
knowledge of the two-body scattering phase shift and 
bound state spectrum (when present). This was sub- 
sequently generalized in the seminal paper of Dashen, 
Ma, and Bernstein [3] where a formulation of statistical 
mechanics in terms of the scattering S'-matrix is given. 
The formulation gives a prescription for calculating virial 
coefficients at any order, and shortly thereafter the be- 
havior of the third order coefficient at low temperature 
was obtained via the S- matrix method [8] . The S*-matrix 
approach to virial coefficients is still an actively research 
topic [9| with recent applications in the field of cold 
atomic gases [To| . However, since determining the S- 
matrix in a general system with multiple particles is a 
high non-trivial task, it is valuable to pursue alternative 
ways of approaching the virial expansion. 

Approximations or assumptions are unavoidable at 
some point. The traditional strategies have been either to 
limit the Hilbert space allowed for the variational many- 
body wave functions or to design schematic Hamiltonians 
aiming for specific features. The latter approach requires 
great care and physical intuition to retain the necessary 
features of the Hamiltonian that will accurately describe 
the phenomena under study. An approach along this sec- 



ond line of reasoning is the use of harmonic Hamiltonians 
where interaction terms are replaced by harmonic oscilla- 
tors. This is extremely convenient from a computational 
point of view as many aspects become analytically ad- 
dressable for both fermionic and bosonic systems [llTflij . 

The replacement of one- and two-body terms in the 
Hamiltonian by harmonic forms leaves the problem of de- 
termining the parameters of the harmonic Hamiltonian 
according to given criteria. Here one must again take 
guidance from physical properties and aspects of the sys- 
tem that are crucial for the system under study. Re- 
cently, we formulated an approach that explores the N- 
body problem in an external parabolic confining potential 
by fixing the two-body interactions to the properties of 
an exactly solvable problem in the same geometry [l5| . 
The two-body information needed is the energy eigenval- 
ues and structural properties of the wave function such 
as radial averages. 

The example studied in Ref. [l5j was short-range in- 
teracting particles in a harmonic trap for which the two- 
body problem can be exactly solved in the zero-range 
limit [16| , Within the field of cold atomic gases, this 
solution was subsequently confirmed by different exper- 
imental groups (l7j ]. The model studied in Ref. [l6| has 
since been used as a start ing point in both nuclear and 
cold atomic gas physics [lg]. Recently, the harmonic 
approximation with parameters fixed to exact two-body 
properties have been applied to particles that interact via 
dipole forces [H|-[23j and shown to accurately reproduce 
numerical few-body results for moderate-to-strong dipole 
strengths ^MM- 

This harmonic method has also been extended to the 
thermodynamical regime at finite temperature [Til [l2l . 
[28j], where it has been shown using a path-integral for- 
malism that the canonical partition function for given 
particle number can be obtained [IJ EH l29j - l3~0 ] . Here we 
consider an alternative approach that uses exact diago- 
nalization of the Hamiltonian and subsequent calculation 
of the relevant degeneracies in the energy spectrum for a 
given number of particles [28] . This is in contrast to the 
usual approximation using the grand partition function 
where only the average particle number is conserved. The 
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method applies to both identical bosons and fermions as 
well as distinguishable particles and combinations of all 
these possibilities [H| Ho, HH ■ The difficult part is to find 
the degeneracies of the iV-body spectrum for the specified 
symmetries required by quantum statistics. The energies 
themselves are easily found from the harmonic oscillator 
solutions. 

Here we consider the virial expansion within the har- 
monic Hamiltonian approximation for identical fermions. 
This paper is a natural extension of Ref. [28[ where the 
thermodynamics of small to moderate size systems of 
fermions and bosons was considered by direct computa- 
tion of the partition function. This requires a numerically 
efficient determination of the level degeneracies which is, 
however, only possible up to moderate particle numbers 
(N ~ 20) even in a harmonic model. The virial expan- 
sion is usually a rapidly converging series and we there- 
fore expect to be able to compute coefficients within the 
harmonic approach. However, there are some subtleties 
with the convergence of the coefficients that must be care- 
fully handled. Therefore we focus almost exclusively on 
the formal development of the virial expansion within the 
harmonic approximation. 

A motivation for our work is the recent investigation 
of universal thermodynamics within cold atomic gas ex- 
periments [32j where the virial expansion has been suc- 
cesfully applied [33[. However, the expansion is general 
and is applicable to other fermionic systems. While the 
typical condensed-matter and cold atom fermion systems 
have two internal (spin or hypcrfinc spin) components, 
we consider the case of single-component fermions here 
in order to keep the formalism simple while still retaining 
the full quantum statistical properties of a Fermi system. 
Multi-component fermionic and bosonic systems will be 
considered in subsequent studies. 

The purpose of the present paper is to formulate and 
explore the harmonic method and the virial expansion 
to prepare for future applications to systems in both 
cold atomic gases, nuclear physics, as well as condensed- 
matter systems. The paper is organized as follows; We 
describe the ingredients of the method in Sec. [TT] The 
numerical illustrations follow in Sec. IIIII for the lowest 
virial coefficients. In Sec. IIVI we summarize and provide 
an outlook for future directions of interest. 



II. THEORETICAL DESCRIPTION 

We first present general definitions of the crucial in- 
gredients. Then we apply the formulation to the results 
of a system of N particles described by a coupled set 
of harmonic oscillator potentials. The approach to the 
large temperature limit is finally modified to exhibit the 
behavior corresponding to the correct high energy spec- 
trum. 



A. Basic definitions 

The classical virial expansion is an expansion of the 
equation of state of a gas of identical particles, usually in 
powers of the number density p, see e.g. [H-Q: 

p = k B Tp{l + B 2 {T)p + B 3 (T)p 2 + ...), (1) 

where p is the pressure, k B is Boltzmann's constant, T is 
the temperature, and the Bi's are the virial coefficients 
of the expansion. In the classical expansion, they are 
related to the intermolecular or interatomic potentials of 
i interacting particles. The advantage of the expansion 
is that it reveals deviations from ideal gas behavior by 
examining just the few-body aspects of the system. The 
leading term is then the ordinary ideal gas expression 
for a non-interacting system. The second term in the 
classical expansion in three dimensions is given by the 
i?2-coefficient 

B2 = -\J [exp(-/?T/ 12 (r)) - l]d 3 r, (2) 

where V%% is the inter-particle potential depending on the 
relative coordinate r, and /3 = \/k B T. The integral in 
Eq. ([2]) is known as a configuration integral, and from it 
one can see that this coefficient only converges for poten- 
tials which decay faster than 1/r 3 . The convergence in 
two dimensions is correspondingly only achieved when V 
decays faster than 1/r 2 . The general coefficient is 



(all independent i-cluster integrals), 

where V is the volume of the system. The i-cluster in- 
tegrals are all the independent clusters containing the 
i-particles, first presented graphically by [34]. A cluster 
in this context can be visualized graphically by thinking 
of i numbered circles with lines connecting them symbol- 
izing the interaction between those particles. These lines 
can be drawn in many different ways, the only require- 
ment for the i-cluster is that all i atoms or molecules 
must be connected to at least one other member of the 
system. The quantum mechanical version of the cluster 
expansion was developed at around the same time by [H] ■ 
For a system of quantum mechanical particles with 
Fermi or Bose statistics, the expansion is most commonly 
performed in the fugacity, z = exp (/3p) , of the system, 
where p is the chemical potential of the A^-body system. 
The grand canonical partition function, Z, is written as 
an expansion in z 

Z = l + zQ 1 + z 2 Q 2 + ..., (4) 

which translates into an expansion for the grand thermo- 
dynamic potential, Q: 

n = -k B TQ 1 [z + b 2 z 2 + b 3 z 3 + ...]. (5) 
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The first virial coefficients can explicitly be written: 

h = {Q2-QDIQ1 (6) 

63 = (Q 3 ~ Q1Q2 + 0?/3)/Qi (7) 

64 = (Q4-QiQ3 + QlQ2 + Ql/2-Qi/4)/Q 1 , (8) 

where Qtv is the canonical partition function for N par- 
ticles of the proper symmetry, that is 



(9) 



where and -Ej are the degeneracy and energy of 
the j'th state of the A^-body system. Thus Qn can 
be calculated solely from the energy spectrum of the 
TV-body system, and all thermodynamic quantities can 
then be obtained from 57 in Eq. ([SJ to the order de- 
sired. The classical virial expansion, Eq. (TTJ, can be 
recovered from the quantum version introduced above 
by using (N) — — §^|t,v and Q = —pV (see for instance 
Ref. [HI where the relation of Bi and bi is also discussed). 
In the present case we have an external trap and the vol- 
ume must be suitable translated into parameters of the 
trap before making detailed comparisons to experiments 
or other studies |36j |. 

In practical calculations, it is more convenient to 
consider the difference between interacting and non- 
interacting systems. We then consider the differences 
AQ n = Q n — Qn and Ab n = b n — b„ \ where the super- 
script (1) denotes a non- interacting system having the 
same A^-body fugacity z. We can then re- write Eq. §5§ 



fi = - k B TQ 1 [Ab 2 z 2 + Ab 3 z 3 



(10) 



where fiW i s the grand thermodynamic potential of the 
non-interaction system with the same fugacity. The dif- 
ferences of the virial coefficients in Eqs. (JB])-© become 



Ab 2 = AQ2/Q1 

Ab 3 = AQ3/Q1 - AQ 2 

Ab 4 = AQ4/Q1- AQ 3 + Q 1 AQ 2 

+ (Qi - (Q 2 1) ) 2 )/(2Q 1 ) . 



(11) 

(12) 
(13) 



These definitions and expressions are general and now 
applicable to a specified set of potentials producing the 
partition functions Q n and the corresponding (differences 
of) virial coefficients. 



B. Harmonic approximation 

The one- and two-body interactions for the TV-body 
system are approximated by second order polynomials 
in Cartesian coordinates. The Hamiltonian is a sum of 
similar terms from each spatial dimension, H = H x + 
H y + H z . In the present work we consider the two- and 



three-dimensional cases. For identical particles of mass 
m the Hamiltonian of the x-direction is 



H, = 



^ 2 ^ 2 1 2 1 \2 

2^ 7^2 + ^mcj m ^ fa ~ x k ) 



dx 



k=l ~ k 
N 



,k=l 



1 2 \ ^ 2 



N(N - 1) 



V„ 



(14) 



fc=i 



where Xi is the coordinate of particle i, ujq and uj in are 
frequencies of the one and two-body interactions, and V sx 
is a constant adjusting the energy to the desired value. 
The factor 1/8 in the second term comes from the use 
of the reduced mass equal to m/2 and to avoid double 
counting in the sum. 

The frequencies and the constant shift can be chosen 
to reproduce certain properties of a modeled system as 
discussed in [lEj]. How these parameters are chosen does 
not affect the computation of the virial coefficients which 
therefore are obtained as functions of the parameters in 
the Hamiltonian. The method is general and applicable 
as soon as a Hamiltonian of the oscillator form is avail- 
able. It is, however, still useful to illustrate by describing 
the procedure for a specific system. We focus on a system 
of identical, spin-polarized fermions confined in an exter- 
nal trap [15| where the fermions interact via a short-range 
potential. Due to the Pauli principle, the particles cannot 
interact in the spherical s-wave channel, and the lowest 
non-trivial interaction will be odd and of the p-wave kind 
(higher odd partial wave channels will be neglected). In 
the zero-range limit, the model of Busch et al. [l6[ can 
still be solved for p-wave interactions in both three- [37} 
and two-dimensional [38| traps. We adjust the interact- 
ing frequency, uji n , to reproduce some property related to 
the spatial extension of the correct two-body wave func- 
tion (the average square radius in the two-body ground 
state in the trap). The shift, V sx , is then added to make 
sure that the exact two-body ground state energy is re- 
produced by the oscillator potential. Notice that this 
(constant) energy shift does not influence thermodynam- 
ics in any essential way. We therefore ignore it for most 
of our discussion except for some comments near the end 
of Sec. |HU 

The accuracy of the harmonic approximation should 
depend on the degree to which the physical two-body po- 
tential allows a quadratic expansion. Naively, this should 
be the case for potentials that have a sizable attractive 
pocket, which is true for many molecul ar p otentials that 
allow a large number of bound states [39| . An example 
is the Morse potential which has been explored in the 
harmonic approximation in Ref. 14]. As mentioned in 
the introduction, for dipolar particles, the harmonic ap- 
proximation is extremely accurate, even in the regime of 
small dipole moments when suitable adjustment of the 
harmonic frequency is performed [19h21| |. In fact, even 
when the real potential is shallow, the energy can be re- 
produced to a few percent accuracy with a careful choice 
of gaussian wave function as shown in Ref. [2l| . 
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For typical cold atomic gas setups, one has harmoni- 
cally trapped atoms interacting via short-ranged interac- 
tions. In the idealized limit of zero-range interactions, the 
two-body problem is exactly solvable as demonstrated 
by Busch et al. [l6j]. In Ref. [HI, the exact solution 
was used to fit the oscillator parameters that provide the 
input for the harmonic approximation. In the strongly- 
bound limit where a deep two-body bound state occurs, 
this choice of parameters leads to the same scaling of the 
energy with particle number that is observed in varia- 
tional approaches pol. |4lj. Also, when the interactions 
have a diverging two-body scattering length (the unitar- 
ity limit), the two-body wave function becomes similar 
to the non-interacting wave function in the trap [HI and 
we therefore expect the harmonic approximation to be 
very good. These features are clearly seen in the one- 
dimensional case as discussed in Ref. [42| . As discussed 
in Ref. [l6| , the one- and three-dimensional cases are very 
similar. The scalings away from these limiting cases are 
similar but not identical to other approaches. In gen- 
eral, we expect the harmonic approximation to give good 
qualitative results for strong interactions, but do not ex- 
pect perfect quantitative agreement with variation or nu- 
merics. On the quantitative side, we note that for one- 
dimensional systems, the three-body energy can be repro- 
duced to within 10 percent in the strongly-bound limit 
[iH . We thus estimate similar accuracy on the third virial 
coefficient. At this point we leave the question of how to 
adjust the two-body parameters of the Hamiltonian and 
proceed with a general discussion for arbitrary parame- 
ters. 

The solution to Eq. (|T4"|) is found by a coordinate trans- 
formation which splits the Hamiltonian into N indepen- 
dent harmonic oscillators with new coordinates and fre- 
quencies related to the normal modes of the A^-body sys- 
tem. For A~ identical particles two new normal mode 
frequencies are produced, that is the external trap fre- 
quency, wo, corresponding to the center of mass motion, 
and the N — 1 times degenerate frequency, ui r , given by 

u? r = NluI/2 + . (15) 

The degenerate frequencies correspond to different types 
of intrinsic (relative) motion which for two particles is 
simply oscillations in the relative coordinate, but in the 
general case correspond to different normal modes of the 
system. The A-body energy spectrum is 

Ej = E cm + E rel + V s , V s = §N(N - 1)V„, (16) 

E cm = huj (n (i) + y) > n oU) = Ei=i n o,i(j), (17) 
E re i = tkJ r (n r {j) + §{N-l)) + V s , (18) 

n r (j) = E£=i Y,k=i nk,i(j) = EiLi n rAj)i (19) 

where j is the index of the A^-body states. no,j is the 
number of excitation quanta in the center-of-mass degree 
of freedom in the i'th direction, while n^ i is the number 
of quanta in the fc'th of the AT — 1 degenerate uj t modes 
in the i'th direction. The energy state j is thus given 



by specifying the number of excitation quanta in each 
normal mode for all dimensions. Here we have divided 
into center of mass and relative contributions since in 
the equal mass case studied here these can be separated 
completely. 



C. The lowest virial coefficients 

The partition function for one particle, Qi, is the triv- 
ial problem of one particle in an external harmonic trap 
of D dimension. Since it has no other particles to inter- 
act with, the spectrum arises only from center of mass 
motion for one particle. The partition function is in this 
case found from Eqs. © and (TT71) to be 

= / cxp(-e /(2T)) y = ^ 6o 
^ VI -exp(-6 /T)/ 2 D K 2T J ' V ' 

where Oo = hui^/kB- 

The partition function for two particles, Q2, is found 
from Eqs. ©, (HH), (|T7]). and CEH]), and can be factor- 
ized into center of mass and relative contributions. The 
center of mass piece is completely symmetric in all the 
coordinates, so it plays no role in determining the overall 
symmetry of the system. It is just a geometric series in D 
dimensions, and in fact equal to Q\, since the frequency 
is that of the external trap. 

The relative motion corresponds to the difference be- 
tween the two individual coordinates. For fermions this 
motion then must provide the antisymmetry of the wave 
function corresponding to an odd number of relative os- 
cillator quanta, n r (j), in Eq. (|18[) . The energies are 
completely specified by the quanta and the related de- 
generacy is easily counted for the relative motion of two 
particles. To obtain an analytical, closed solution, it is 
convenient to consider the individual Cartesian quanta. 
In 2D, one merely needs to keep either quanta, n r >x or 
rir^y odd in x or in y-directions. In 3D, there are two pos- 
siblities, in that all three quanta are odd or two are even 
and one is odd. With these restrictions the summation 
in Eq. © leads for D = 2 to 

fe = «'(T^OT)) aexp( - 6s/r) ' (21) 

where 8s = VsI^b and 8 r = hjj r /ks- For D — 3, we 
find instead 

Q 2 = Q 1 cxp(-e s /r) (22) 
3 exp(-56 r /(2T)) + exp(-96 r /(2T)) 
X [1 - exp(-2e r /T)]3 

The next terms involve the partition function, Q3, for 
the three-body system. A completely closed form solu- 
tion for the partition function is not found, and we keep 
the expression as an energy sum over three-body states. 
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The center of mass summation is again performed ana- 
lytically, and Eqs. (|16p and (|T5)) lead to the other factors 
amounting in total to 



Q 3 = Qiexp(-e s /T)}^.g ; exp(-(Z 
1=0 



D)Q r /T), (23) 



where the summation over all states is reduced to run 
over all integers. The difficulty is then only to know the 
corresponding degeneracy, gi . This number of states of a 
given excitation energy is found by the method described 
in [28|. The expression in Eq. (|23[) is formally the same 
in two and three dimensions but the degeneracy factors 
differ substantially. The infinite sum must in practice be 
truncated at some level of excitation. Ideally this is after 
convergence is reached. However, this depends strongly 
on the value of the temperature and we therefore first 
must decide how large T-values we need to investigate. 

The differences between interacting and non- 
interacting virial coefficients are found from Eqs. (fTTj) 
and (|12|) . The non-interacting partition functions are 
structurally the same as the interacting ones, only with 
8 r replaced by Oo, and 8 s = 0. Using the expressions 
in Eqs. (|2"T|) . (f2"2"j). and (|2"3")l , we can therefore easily find 
A62 and A63. For A^2 in two dimensions we explicitly 
get 



Ab 2 = cxp(-e s /T) 



exp(-8 r /T) 
1 - cxp(-26 r /T) 



exp(-9 /r) 
1 - exp(-2e /T) 



(24) 



and a slightly more complicated expression for three di- 
mensions from Eq. (1221) . In the same way we can of course 
re- write A63 in terms of the expression for Q3 in Eq. (1231) , 
but due to the lack of closed form it does not provide any 
further information. 

Notice that the low-temperature limit of the virial co- 
efficients above is determined by the value of the shift, 
6s- From Eq. JM]) we see that Ab 2 will vanish for T -> 
whenever we have 8s + 28 r > 0, similarly for the higher 
virial coefficients. In the following we will mostly dis- 
cuss the case 85 = 0, i.e. no shift at all, since this is 
the relevant case for thermodynamics. For generality, we 
comment briefly on the influence of the shift for all tem- 
peratures in Sec. IIIII 



D. Large-temperature limits 

At high temperature, the interacting system should ap- 
proach a non-interacting system as the kinetic energy 
dominates the potential. This does not imply that all 
Abi must vanish at large T, since deviations in the low- 
energy spectrum between interacting and non-interacting 
systems easily produce different large-temperature con- 
tributions. This can be seen by dividing the sum over 
states in the partition functions in low- and high-energy 



parts. Even if we assume that the high-energy interact- 
ing and non-interacting spectra become equal and would 
thus not contribute to Abi, the low-energy parts remains 
different and this yields a contribution to the virial coeffi- 
cient for all temperatures. However, such difference must 
remain finite since it arises from a finite energy interval. 

As stated before, the virial expansion does not work 
for potentials that do not vanish at a fast enough rate, 
so one would think that the harmonic potential, which 
does not vanish at all, would cause problems. Indeed, if 
we use the derived expressions all our Abi diverge with 
increasing temperature. The origin of this problem is 
simply that the energy spectrum is obtained for the N- 
body solution of a temperature independent Hamiltonian 
adjusted to reproduce ground state properties. This does 
not account for the influence of temperature on the effec- 
tive interactions and, in turn, on the energy spectrum. 
This may also be expressed in terms of an excitation en- 
ergy dependence of the effective interaction as seen for 
example in the variation of the mean-free-path. Since 
excitation energy on average can be related to tempera- 
ture these formulations are equivalent. 

To pin-point the problem and subsequently cure it we 
start with A&2 ■ In the limit of high temperature, Eq. ([24]) 
can be expanded to leading orders in T to give 



Ab 2 {T -> 00) = T 2 



1 



1 



462 48 2 



T 



492 



Equivalently we find for three dimensions that 



Ab 2 (T ^00) = T 



1 



1 



(287 2(-) ( ; 



8g 
283 



(25) 



(26) 



T/e| 15 



15 

e~ 



Obviously Ab 2 diverges as T to the power of the dimen- 
sion D, unless of course interacting and non-interacting 
frequencies are precisely equal. The next orders are in- 
dependent of T or vanishing with increasing T. 

The unphysical divergence is here clearly seen to orig- 
inate from the difference between the energy spacing 
found by a fit to ground state properties and the spac- 
ing for a non-interacting system. The proper thermo- 
dynamics requires a correct spectrum for all energies or 
temperatures. 

In order for A62 to be finite at large temperatures, 
io r must approach ujq. From Eq. (|15[) this seems most 
reasonably achieved by a vanishing interaction frequency 
ujin. Furthermore, Eqs. (|2"5"|) and (|2"B1) imply that the 
energy shift adjusted to fit the ground state energy also 
must vanish in the large-temperature limit. To get finite 
A62 we introduce cutoff functions F(T) and G(T) into 
Lo in and V s : 



N 



>jj,. 



V s -> V S G{T) 



(27) 



Without further constraints, there is a great deal of free- 
dom in the form of F(T) and G(T), but they must satisfy 



6 



the conditions of approaching zero respectively as 1 /T D 
and at high temperatures. The simplest such 

functions 

(D-l) 



F(T) 



T 



T + T 



> G(T) 



T 



T + T 



, (28) 



where T is a cut-off parameter which indicates at what 
temperature the spectrum continuously should shift to 
the non-interacting spectrum. These functions regularize 
the high temperature behavior and now A62 approaches a 
constant at high temperatures. The limits in the different 
dimensions can be found from Eqs. (|25[) and (|26[) by use 
of Eqs. ([27]) and ([28]). that is 



Ab 2 (T -> 00) 



Nu? n (k B T ) 2 Vsk B T 
' 8 (Huj ) 2 4(/iw ) 2 



(29) 



for two dimensions and for three dimensions we get 



Ab 2 (T -> 00) 



3N cof n (k B T ) 3 V s {k B ToY 



8 w 2 (nw ) 3 2(hu} ) 3 



(30) 



where Af = 2 for this second virial coefficient. The high 
temperature limit is then determined by the interaction 
frequency, ground state energy shift, and the cut-off tem- 
perature. The two terms have opposite sign since Vs is 
negative, and initially introduced to balance the zero- 
point energy of the oscillator such that large uji n is cor- 
related with a large negative Vs- 

Once the A&2 is regularized we turn to A&3 which 
is more complex as seen in Eq. (TT2]) . It consists of 
two terms, the difference between interacting and non- 
interacting partition functions for three particles and 
AQ2 = QiA&2- This latter term contains the already 
regularized factor, A6 2 , but Qi in Eq. ([20]) also di- 
verges as — > T d /Qq at high temperature. Therefore 
A63 can only remain finite in the limit when the differ- 
ence (Q3 — Q^)/Qi diverges precisely as Q\. 

However, this is not even sufficient because the T D di- 
vergence from Q\ leads to divergence of the terms in A6 2 
vanishing as 1/T in two dimensions and as 1/T and 1/T 2 
in three dimensions. These terms all must be cancelled by 
corresponding diverging terms in (Q 3 — Q3 )/Qi- These 
conditions seems impossible to meet, but nevertheless 
this miracle seems to occur. We find that A&3 also is 
regularized with precisely the same cut-off function as 
used for A&2 and this occurs in both two and three di- 
mensions. Problems with divergences in A&3 have been 
reported by other authors (43| . They resolved it by sep- 
arating the system into a 2+1 subsystem in which the 
divergence was then removed. While we do not need to 
do this to obtain a finite A63, the underlying deep reason 
for this is still unclear to us. 

The form of the cut-off in Eq. ([2"7]) can be changed in 
many ways. The transition can set in at different tem- 
peratures and be more or less fast. We have chosen to 
use only one parameter, To, but we tested with another 
functional form: 



which has the same high temperature behavior. Both 
A6 2 and A63 are again regularized in the large- 
temperature limit. The values of A6 2 and Ab 3 are both 
larger in magnitude at all temperatures simply because 
this cut-off function is larger for all T. 

Instead of the temperature cutoff in Eq. (|2"T]) , we could 
choose a cutoff in excitation energy. This might at first 
appeal as being more physically reasonable since this di- 
rectly amounts to changing the A-body spectrum at high 
excitation energy to the non-interacting iV-body spec- 
trum. However, this conclusion is rather shaky. Com- 
pletely different configurations specified by sets of quan- 
tum numbers can give precisely the same energy. This is 
easily seen in a single particle picture by comparing (i) 
states of the same total energy arising from a few parti- 
cles at very high-lying levels and all others in the lowest 
possible levels with (ii) the opposite where all particles 
are in levels of intermediate energy. 

Thus, a given excitation energy already corresponds to 
an average over many configurations very similar to the 
configurations occupied for a given temperature. In any 
case we investigated a cutoff in excitation energy instead 
of temperature, that is expressions of similar simplicity 
as those in Eq. 



F(E) = 



En 



E + E 



G(E) = 



En 



E + E 



(13-1) 



, (32) 



F(T) = (1 - exp(-T /T)) 



D 



(31) 



where E is the excitation energy. This function is in- 
serted in Eq. ([27]) in the place of F(T), and when the 
partition function is being calculated as the sum over 
states, the step in energy taken between states is now 
dependent on the position in the spectrum. The mean 
field spacing is reached for energies higher than Eq- Im- 
plementing this cutoff successfully regularizes A62, but 
to remove the divergence for A6 3 , it is necessary to use 
a constant, Eq, different from that of Ab 2 - No obvious 
relation is found and generalizations to higher Abi seems 
to be at least very impractical. 



III. NUMERICAL RESULTS 

The virial coefficients can now be calculated numeri- 
cally. In the model they are completely determined from 
trap frequency U)q, interaction frequency u>i n , shift en- 
ergy Vs , and cut-off function and related parameters To . 
We use the trap frequency as energy unit which implies 
that results for any other value of ujq can be obtained 
by scaling all energies k B To, and Vs by fiwo- The 

energy shift, Vs, is introduced to adjust to the correct en- 
ergy and has no effect on eigenvalues and corresponding 
wave functions. It is strongly dependent on which model 
is approximated. We shall therefore first investigate the 
general dependencies on u>i n and To which in turn can be 
related to specific models. Afterwards we shall separately 
investigate the dependence on the shift. 
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FIG. 1: A&2 in 2D as a function of temperature with a model 
space of 160 oscillator shells and zero energy shift. The cutoff 
used was the function in Eq. (|28[) . The values of To and u)i„ 
(in units of wo) from bottom to top are (To,uj in ) — (0.51, 5.0), 
(0.25,5.0), (0.13,5.0), (0.25,2.5), (0.13,5.0), and (0.062,2.5). 
The To values are cji„/27r 2 and twice and half that value. 



< 




5 10 

k B T/nco 



FIG. 2: Same as in Fig. [T]for A63 in 2D. The parameters are 
also the same but must be read from top to bottom. 



A. Virial coefficients 

The cut-off parameter is essential for the behavior of 
the expansion coefficients. The size of an appropriate 
value can be estimated by inspection of the effect it is 
designed to simulate. The first excited state appears at 
an excitation energy of huj r which is a single quasi par- 
ticle excitation. Therefore Huj r represents a shell gap 
to be overcome by thermal excitations. This gap is 
known to wash out at a critical temperature, T given 
by fcsT27r 2 w fiw r , see [44|,|45[. This surprisingly large 
factor, 2n 2 w 19.7, suggests a rather small relative value 
of To proportional to u> r , that is ksTg » huj r / (2ir 2 ). 

To illustrate the dependence we show A&2 in Fig. [T] as 



function of T for different interactions and cut-off values. 
We see the general behavior of a second order increase 
from zero at T — 0, and the smooth curvature before 
bending over to reach the saturation value. The expan- 
sion coefficient is a rather strongly increasing function of 
both interaction frequency and cut-off parameter. The 
functional form of the cut-off function in Eq. ([28)) implies 
that the saturation value as well as saturation tempera- 
ture both depend rather strongly on To. 

The overall behavior must be understood in the model 
even at uninterestingly large temperatures. We continue 
to show A63 in Fig. [2] for for two dimensions for the same 
set of parameters as in Fig. [TJ Qualitatively the same be- 
havior except for the overall opposite sign. However, the 
temperature dependence is faster, the saturation values 
are larger, as seen for the small interaction frequency with 
the small To . For higher values we observe a tendency to 
form a fiat region which quickly becomes an increasing 
function. 

At low temperature, the coefficients vanish since both 
the interacting and non-interacting partition functions go 
to unity at low temperature (in the abscence of any en- 
ergy shift). This might seem odd since our A&^'s are the 
difference between a non-interacting and interacting sys- 
tem, which should increase at low temperatures when the 
interactions are more significant compared to the kinetic 
energy. Since both of the partition functions are small 
at low temperature, the only signature we see of this is 
that the relative differences of the partition functions, for 



example (Q2 — Q^) IQ^' > does increase. 



(i) 



.0 

< 



< 




FIG. 3: Virial coefficients A&2 < and A&3 > as a function 
of interaction frequency uii„ for three different temperatures. 
Notice that from Eq. (|15[) we always have un n > \/2luq. On 
this plot, ksTo/hwo = 0.25 for all curves except the short- 
dotted line which has ksT/huj — 5.0 and fesTo = fcu}i n /2n 2 
(adjusted at end point on the horizontal axis). 

The effect of the interaction frequency is plainly to in- 
crease the virial coefficient, which could be seen in Figs.Q] 
and [2 A more precise dependence on interaction fre- 
quency can be seen in Fig. [3J The coefficients vanish 
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FIG. 4: A63 in 2D as a function of temperature for three 
different model space sizes. The cut-off function is that given 
in Eq. (|28[l . The other paramters are fcsTb = 0.25fiwo and 
(Jin = 2.5wo 

for small uii n , as then there is no difference between the 
non-interacting and interacting system, but then increase 
rapidly, especially after <jj„ > loq. The increase is more 
dramatic at higher temperatures, which are closer to the 
saturation value. It does appear, however, that for any 
temperature the behavior of the coefficients is faster than 
linear. 



B. Hilbert space and cut-off function 

The apparent lack of saturation for A63 at high tem- 
peratures is a very unsatisfactory feature. Fortunately, 
it seems to be an effect of the model space truncation at 
high excitation energies in the calculation of the partition 
function. This happens because of the subtle nature of 
the cancellation which removes the divergence and leads 
to saturation. The piece, Q\, in AQ 2 in Eq. (fT2j) contains 
all states to infinitely high excitation energies since it is 
calculated analytically. This piece eventually overwhelms 
the AQ3 term in Eq. (|T2j) which is calculated numer- 
ically and consequently arises from a truncated energy 
spectrum. 

We demonstrate this in Fig. @] where we show A63 for 
three different model space sizes. The higher energies 
we include in the numerical calculation, the larger is the 
region of the flat saturation interval, and the larger tem- 
peratures before the divergence sets in. This is very reas- 
suring allowing us to ignore the unphysical region of all 
temperatures above the flat region. 

One uncertainty in the method to recover the high- 
energy non-interaction limit is the function describing the 
disappearance of shell effects. In Fig. [5] we illustrate this 
dependence of the virial coefficients by results from use of 
different cutoff functions, that is the rational expression, 
Eq. ([25)) , and the exponential functon, Eq. (|5T|) . The 



FIG. 5: Virial coeffecients in 2D as a function of temperature 
for different cut-off functions. The rational function cutoffs, 
Eq. (|28[) . are labeled with an 'R', and the exponential cutoff, 
Eq. (|31|) . with an 'E'. uii n = 2.5o;o was used for all of the 
curves. 



virial coefficients using the exponential cutoff are larger 
for all temperatures before finally merging into the same 
high temperature limit. This is due to the larger cut-off 
function at all temperatures which leaves the reduction 
to take place somehwat faster at larger temperatures al- 
though the same final saturation is reached for T much 
larger than Tq. 
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FIG. 6: A&2 (lower solid line) and A&3 (upper dashed line) 
in 2D as a function of temperature examining the effect of 
the energy shift. The coeffecients behave similar to what was 
shown previously except at very small temperature. For this 
plot, LJi„ = 2.0o;o with V s = — 5.8Ifiwo, and fcsTb = 0.5/kj. 
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C. Effect of the energy shift 

The energy shift in the Hamiltonian has a different ef- 
fect on the virial coefficient. We illustrate by the exam- 
ples in Fig. [5] The large temperature limit remains finite 
by construction as shown explicitly in energies Eq. (|29|) 
for one case. Otherwise the behavior at large T is very 
similar to that of zero shift, where convergence is essen- 
tial or at least a flat region at high T is necessary. The 
energies enter in the partition function in the exponent. 
Contributions disappear from energies much higher than 
the temperature. 

However, at very small temperatures the same con- 
tribution from the shift can produce unphysical results. 
This is seen at very low T in Fig. [5] where a narrow and 
large peak is present in both virial coefficients. The nu- 
merical reason is obvious since the negative shift divided 
by a small temperature produce a very large value. This 
occurs for temperatures much smaller than ujq and much 
before the statistical treatment is meaningful. 

The negative shift is not necessarily sufficient for a di- 
vergence at low temperatures. The shift energy must 
completely eliminate the zero-point energy, and thus for 
polarized fermions Vs < —DNftw r . The parameters in 
Fig. [S] provide a shift large enough in magnitude to give 
a spike towards large positive values at very low tem- 
peratures. This changes the sign of A62, as the term 
containing the shift eventually (small T) becomes larger 
than the non- interacting term in Eq. 

For A63, the singularity at zero temperature comes 
at a slightly higher temperature since the shift energy 
is three times larger than for two particles. Again this 
occurs for temperatures smaller than those allowing a 
statistical treatment. 



D. Three dimensions 

We have so far only shown results for two dimensions. 
The method is, however, as applicable in three dimen- 
sions where the coefficients look qualitatively similar to 
those in two dimensions, see Fig. [7] We notice the in- 
crease from zero at T = 0, then a decreasing derivative 
resulting in saturation or in a tendency towards satura- 
tion, and finally the divergence at large temperature for 
A63 due to mismatch between the numerical Q3 and the 
analytical Q\. The results are more sensitive to the cut- 
off parameter because higher powers enter in two than 
in three dimensions, as seen for example by comparing 
Eqs. (H!]) and §U§. Indeed, if one uses T = u} m /(2n 2 ), 
then the magnitude of the coefficients is quite small, of 
order 10~ 2 . Also, the truncation effect in A63 is more 
sensitive to the cutoff parameter and can begin at com- 
paratively low temperatures when compared to 2D, see 
Fig.H 

The virial expansion is most efficient when the size 
of the coefficients decrease with the order. Therefore, a 
requirement of A63 < A62, leads to a condition on the 




-1 - 



5 10 15 

k B T/neog 

FIG. 7: The virial coefficients in three dimensions, A62 < 
(lower lines) and A63 > (upper lines), as function of temper- 
ature. The cut-off temperatures, To, and u)i n is varied among 
the different curves. On the A&2 side (below zero) from bot- 
tom to top it is (fcsTo/fttiJo, uJin/uo = (1.00,2.5), (0.50,3.5), 
(0.50,2.5), (0.50,1.5), and (0.25,2.5). For A63 (above zero), 
the lines have the parameters just listed but in the inverse 
ordering. For this plot, there is no energy shift Vs — and a 
model space of 160 shells is used. 



maximum size of To. This demand is more restrictive for 
three than for two dimensions. 



IV. SUMMARY AND OUTLOOK 

We have discussed the virial expansion technique and 
a quantum mechanical formulation is sketched from an 
analogous classical expansion. We apply the formu- 
lated method to a harmonic approximation to the N- 
body problem for identical fermions. A key step in 
this approach is the adjustment of the harmonic one- 
and two-body parameters to pertinent properties of the 
corresponding two-body problem that holds information 
about the exact interaction which is approximated by 
a harmonic form. Once these are obtained, the result- 
ing A-body Schrodinger equation may be solved exactly 
and the spectrum can be used to compute the partition 
function. The second and third order virial expansion 
coefficients are obtained by direct calculation of the two 
and three-body partition functions. Here we are inter- 
ested in the details of the formal development of a virial 
expansion and we therefore vary the parameters of the 
harmonic interaction terms freely to study the behavior. 
The mapping to realistic two-body properties is straight- 
forward. 

The virial expansion can be reformulated in terms of 
deviations between non-interacting and interacting sys- 
tems. Importantly, the virial coefficients have an un- 
physical divergence for large temperatures. It arises in 
the formulation because the increasing temperature pop- 
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ulates higher and higher excited states. Their average 
properties can be very far from the ground state prop- 
erties, and eventually the results should resemble those 
of the non-interacting system where the kinetic energy is 
decisive. This is achieved by modifying the energy spec- 
trum by a function of temperature smoothly connecting 
the low-temperature ground state dominated and high- 
temperature non-interacting spectra. 

To achieve the goal of removing the divergence, the 
modification function must reduce the initial interaction 
frequency to zero by a large-temperature behavior where 
the power of the temperature is equal to the spatial di- 
mension of the system. The divergence is removed from 
the second order expansion coefficient and it turns out 
that the same modification function removes the diver- 
gence from the third order term. This result is highly 
non-trivial since the third order divergence is of a very 
different origin from that of second order. This is em- 
phasized by an attempt to use an energy (in contrast 
to temperature) dependent modification function which 
removes the second order divergence but requires an ad- 
ditional adjustment to remove the third order divergence. 
The temperature modification is then the more promis- 
ing approach for applications where higher orders have 
to be calculated. 

We find that the critical temperature value describing 
where the adjustment of the spectrum should take place 
and the rate of the modification should be about 2tt 2 
smaller than the two-body interaction frequency. This is 
an analogy to the smearing of shell effects by temperature 
in an A^-body finite system described by its single-particle 
spectrum. This function is exponential and the rate is 
precisely the single- particle frequency divided by 2ir 2 . In 
our case the modification function is chosen to have a 



rate of change of the iV-body energy spectrum which is a 
rational function of temperature to a power that depends 
on dimension. However, the modification takes place on 
the energies appearing in the exponent of the partition 
function. The precise shape of the temperature modifi- 
cation function is not essential for the overall properties 
unless one consideres extreme cases. A sensible choice 
of modification function that regularizes the divergence 
will thus yield a formalism that can be adjusted to low- 
energy properties and make further predictions for the 
many-body problem. 

While we have focused on two- and three-dimensional 
systems in the current presentation, a good testing 
ground for the formalism discussed would be some of 
the exactly solvable models that are known in one- 
dimensional systems 46|. A good example is the N- 
boson problem with zero-range interactions studied by 
Lieb, Linniger, and MacGuire [l?) for which the har- 
monic approximation can be applied 42], or to dipolar 
molecules in one-dimensional setups where Af-body clus- 
terized bound states can easily form |48j |. In the strong- 
coupling domain these should be well-described within a 
harmonic approximation approach. 

In summary, we have demonstrated that the harmonic 
approximation employed at the Hamiltonian level gives 
a divergent set of virial expansion coefficients that must 
be regularized at high temperature. This can be done by 
a careful choice of temperature-dependent spectral mod- 
ification that will render the virial coefficients finite at 
all temperatures. The ease of solving the harmonic N- 
body problem and subsequently calculating the virial co- 
efficients makes this an attractive approach to compute 
many-body properties. 
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